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ABSTRACT 

The propagation trajectories of the highest energy cosmic rays (HECRs) are deflected by not only intergalactic 
magnetic field but also Galactic magnetic field (GMF). These magnetic fields can weaken the positive correlation 
between the arrival directions of HECRs and the positions of their sources. In order to explore the effect of 
GMF on the expected correlation, we simulate the arrival distribution of protons with the energy above 6 x 10'^ 
eV taking several GMF models into account, and then test the correlation between the protons and their sources 
assumed in the simulation. The dependence of the correlation signals on GMF models are also investigated. The 
correlation can be observed by accumulating ^ 200 protons in a half hemisphere. Typical angular scale at which 
the positive signal of the correlation is maximized depends on the spiral component of GMF models. That angular 
scale is ^ 5° for bisymmetric spiral (BS) GMF models and ^ 7° for axisymmetric spiral (AS) GMF models if 
the number density of HECR sources, n,, is ^ 10""* Mpc"^. An additional vertical (dipole) component of GMF 
affects these angular scale by 0.5° - 1°. The difference between the correlation signal for the BS models and that 
for the AS models is prominent in the northern sky. Significance of the positive correlation depends on source 
distribution. The probability that the number of simulated HECR events correlating with sources is smaller than 
the number of random events correlating with the same sources by chance is much less than 10"^ 3(j) in almost 
all the source distributions with = 10"^ Mpc"-* under 200 protons detection, but ~ 10% of source distributions 
predicts the chance probability more than 10"^ in the AS GMF model. In addition, we also briefly discuss the 
effect of GMF for heavy-nuclei dominated composition. 

Subject headings: cosmic rays — methods: numerical — ISM: magnetic fields — Galaxy: structure 



1. INTRODUCTION 

The origin of the highest energy cosmic rays (HECRs) is 
one of the biggest mysteries in astrophysics. Although theo- 
retical candidate s of HECR s ources have been proposed (see 
iBhattachariee & Sigll (l2000t) : iTorres & Anchordoqui (2004) 
for reviews), HECR sources have not been successfully iden- 
tified. A difficulty to identify HECR sources is the existence of 
Galactic magnetic field (GMF) and intergalactic magnetic field 
(IGMF). Since HECRs are charged, these magnetic fields de- 
flect the trajectories of HECRs from their sources to the Earth, 
i.e., the arrival directions of HECRs are shifted from the posi- 
tions of their sources. 

Recent progress to construct HECR observatories with larger 
effective area has enabled us to accumulate a number of 
HECR even ts and to draw aniso t ropy in the arrival distributio n 
of HECRs (iTakeda etal. l ll999t [Abraham et ani2007l l2008h . 
The anisotropy includes information on the positions of their 
sources and intervening magnetic fields from the sources to 
the Earth. Also, the Pierre Auger Observatory (PAO) found 
the correlation between the arrival directions of 27 HECR 
events above 5.7 x 10'^ eV and the positions of extragalac- 
tic astrophysical objects in local Universe within z = 0.018 
with angular scale of 3.1° (Abraham et al. 2007, 2008 ). Sev- 
eral au thors have also analyzed the PAO events. Georg e et al. I 
(l2008h showed significant correlation of the PAO events with 
hard X-ray selected active galactic nuclei (AGNs) by us- 
ing the 2 dimensional generalization of Kolmogorov-Smilnov 
test. Generally, the correlation between the PAO events and 
galaxies, which are a representative of matter distribution 
in local Universe, has been shown (iGhisellini et al. I l2008t 
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iKash ti & Waxman Il2008t iTakami et al. Il2009h . iTakami et al.l 
(2009) estimated the typical angular scale of positive corre- 
lation between HECRs and galaxy distribution and bounded 
the deflection angles of HECRs in local Universe. All these 
works indicate that the origin of HECRs is extragalactic ob- 
jects distributed with being comparable to galaxies. Further- 
more, in an optimistic view, we can interpret these results 
as the evidence that intervening magnetic fields are not so 
strong that the arrival directions of HECRs lose information 
on the posi tions of their sources as predict ed by Dolag e t al. I 



('2005') and 'Takami et al. . , _ 

Kotera & Lemoine ( 20081) iDas et al. I (|2008|) and Ryu et aL 



(120061) (but see ISigl e t al 



(|2009) for strong IGMF). This interpretation implies that a 
more number of HECR events can unveil the positions of their 
sources as strong event clusters with i n sma ll angular scale, as 
demonstrated by Blasi & de MarcoJ ( |2004|) without magnetic 
field. Takami & Sato (2008a) showed that such a signal is also 
predicted for the highest energy protons (HEPs) even consid- 
ering a structured IGMF model, which reflects the matter dis- 
tribution of local Universe actually observed. We can regard 
these works as predictions of the PAO result. However, it is still 
under debate that astrophysical objects which correlate with the 
PAO events are really HECR sources because these objects have 
weaker activity than theoretically predicted source candidates 
dMoskalenko et al. 2009). 

However, a recent PAO re sult showed tha t the si gnificance 
of the correlation reported in [Abraham et al. I (l2007l 12008.) de- 
cre ased, though the cor relation is still significant at more than 
2a (iHague et al. 11200 9*). In addition, the High Resolution Fly's 
Eye (HiRes) reported that their is no correlation between HiRes 
events and extragalactic objects based on a similar analysis to 
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the analysis of the PAO dAbbasi et al. II2008I) . Thus, we should 
reconsider the correlation between HECRs and their sources, 
and therefore the effect of intervening magnetic fields. The 
large deflection angles of HECRs could weaken the correla- 
tion. The possibilities that the deflection of HECR trajectories 
is larger than one expected are (i) GMF affects the deflections 
significantly (ii) the main composition of the HECRs is heav- 
ier than protons (iii) the IGMF model used in Takami & Sato I 
(l2008a ) was optimistic, i.e., the deflection of HEPs in this 
IGMF model is smaller than that in the real Universe. In this 
study, we focus on the possibility (i) assuming that the compo- 
sition of HECRs are purely protons and IGMF subdominantly 
affects the deflections of HEPs. The other possibilities will be 
discussed in ^ 

Cosmic rays arriving at the Earth are inevitably af- 
fected by GMF. This importance has driven us to con- 
sider the gropagation of HECRs in t he Galac tic space 



; pro p 
199^ 
et al. 



(IStanev II1997I: [M edina-Ta nco et al. Ill998t iHarari et aL Hl999l: 



2001; Alvarez- Muniz et al. 120021: iHarari et al. I 



Kalashev i 

2002allbt iTi nvakov & Tkachev I 120021: iProuzaetal. 



^shiguchi et al. 2003b, 2004; Tinvakov & Tkachev 
Kachelriess et al. 2007; Takami & Sato 2008b; Cuoco et alj 
2008: Golup et al. 20091). Many of these works have mainly in- 
vestigated the propagation of HECRs itself in detail. Based on 
the results of these studies, the deflection angles of the trajec- 
tories of HEPs are highly dependent on their arrival directions 
reflecting the structure of GMF and can be less than 10° at 
the highest energy (~ 10^'' eV) except in the direction of the 
Galactic center. The angular dependence of the deflection of 
HEPs distorts the apertures of HECR detec tors to extragalac- 



tic spa ce, so-called magnetic lensin g effect dHarari et al. II 19991 
l2002allbl: iKachelriess et al. 1 12007|) . Some of them guessed 
HECR sources using experim ental data and HECR trajectories 
calculated in GMF models. 'Golu p et alH ( |2009|) studied the 
possibility to reconstruct the positions of HECR sources and 
several pr operties of intervening magn etic fi eld from detected 
HECRs . lAlvarez-Muniz et al. I ( l2002h and lYoshiguchi et al. I 
(I2003bh calculated the arrival distribution of HECRs assuming 
source models, and then analyzed the arrival distribution. How- 
ever, they mainly investigated the auto-correlation of HECRs. 
The cross-correlation between HECRs and their sources was 
hardly considered. 

In this study, we investigate the cross-correlation between 
the arrival directions of HEPs and the positions of their sources 
containing the effect of GMF. For this purpose, we simulate the 
arrival distribution of protons above 6 x 10'^ eV taking their 
propagation in magnetized Galactic space into account. Since 
the positions of sources of simulated HEPs are known in sim- 
ulations, we can calculate cross-correlation functions between 
HECRs and their sources. The modelling of GMF is still con- 
troversial, as different models of GMF have been adopted in the 
previous studies introduced above. Thus, we also examine the 
dependence of predicted correlation signals on GMF models. 

This paper is laid out as follows: ^is devoted to explain 
GMF models, a HEP source model, the calculation method of 
the arrival distribution of HEPs, and a statistical quantity to use 
our correlation analysis. |3]is dedicated to giving calculation 
results. First of all, we investigate the effect of GMF to the 
propagation trajectories of HEPs in ^3.1! Then, we calculate 
the cross-correlation between the arrival directions of simulated 
HEPs and their nearby sources, and examine the angular scale 
at which positive correlation is predicted in 33.21 In 33.31 we 



estimate the angular distance within which positive correlation 
signals are maximized and its significance. Finally, we summa- 
rize results obtained in ^and briefly discuss the possible effect 
of IGMF and heavy-nuclei dominated composition. 

2. CALCULATION TOOLS 

2.1. GMF Models 

GMF can be divided into 3 components; a spiral component 
following the optical spiral arm of the Galaxy, a vertical com- 
ponent, and a turbulent component. These components of GMF 
are modelled as follows. 

The spiral component could be well described by axisym- 
metric (AS) or bisymmetric (BS) field models. BS models have 
several field reversals, while AS models do not. Although a 
field reversal at ^ 0.5 kpc inside from the solar system is well 
established, the existence of the other field reversals predicted 
by BS models, i.e., which is a better modelling, is still contro- 
versial. Recognizing this uncertainty, we adopt both BS and 
AS models in this study and discuss the dependence of the po- 
sitional correlation between HECRs and their sources on GMF 
models. We u se the parametrization of both GMF models by 
IStanevI (Il997l) . described below. 

The radial and azimuthal components of a spiral magnetic 
field in the Galactic plane are given by 

=B(r||,6'||)sin/?, Be = B{r\\,0\\)cosp, (1) 
where r|| and 9\\ are the galactocentric distance and azimuthal 
angle around the Galactic center, respectively. Note that ^?|| is 
defined as increasing clockwise and 6'|| =0° corresponds to the 
direction of the Galactic center, p is the pitch angle of the spiral 
field in the neighborhood of the solar system, set to p = -10°. 
The field strength at a point {r^^,d\\) in the Galactic plane is 
written as 

f /7(r||)cosf6'||-;91n^') : BS 

^('■Ih^l|)= V (2) 

1 b{r\\) cos(6i||-/31n^j : AS. 

Here /3 = (tan = -5.67 and ro = 10.55 kpc is the galactocen- 
tric distance of the location with maximum field strength at 0| | = 
0°, which can be expressed as ro = (/^^ + c/)exp [-(7r/2)tan/?] , 
where = 8.5 kpc is the distance of the solar system from the 
Galactic center and d = -0.5 kpc is the distance to the nearest 
field reversal from the solar system in the BS model. Negative d 
means that the nearest field reversal takes place in the direction 
of the Galactic center. Note that the AS model does not have 
any field reversals. b(r\\) is the radial profile of the magnetic 
field strength modelled by 

Mr||) = B„^, (3) 

where Bq = 4.4/iG, which corresponds to 1.5/iG in the neigh- 
borhood of the solar system. In the region around the Galactic 
center (r| | < 4 kpc), the field is highly uncertain and therefore 
assumed to be constant and equal to its value at r|| =4 kpc. The 
spiral field is assumed to be zero for r|| > 20 kpc. 

For the spiral halo field, we adopt an exponential decrea sing 
model with two scale heights varez-Muniz et al. l2002h . 
m n I h J ^'^P^-I^l^ " W< 0-5 kpc 



) : \z\ > 0.5 kpc 



(4) 

where the factor exp(-3/8) makes the field continuous on z- 
The parity is represented as 

B(r\\,9\\,z) : S type parity 
-'B('"||>^'||)^) ■ A type parity. 



5(r||,^||,-Z): 
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The combination between the 2 spiral structures and the parity 
produces 4 different spiral GMF models. 

There is little observational evidence of the vertical com- 
ponent of GMF. Near the solar system, a vertical component 
of magnetic field with strength of 0.2 - 0.3 ^(G, which is di- 
rected toward the north Galactic pole, is observed (i Han & Qiao I 
119941) . Also, a lot of nonthermal gaseous filaments perpendic- 
ular to the Galactic plane with tens of /iG to mG have been 
discovered near the Galactic center (.Han 2007) . Little infor- 
mation allows us to construct the re alistic model of the vertical 
component of GMF. ^ Stanev I d 1 997h assumed uniform magnetic 
field with z-direction which is a possible effect of the existence 
of Galactic wind. Here, we a dopt a dipole field as a vert ical 
component of GMF following iAlvarez-Muniz et al. I ( l2002h . A 
dipole field is described as 

=-3/iG sin(/)cos0cos(/3/r-^, 
By = -3/XGsin(/)COS(/)siniy9/r-', (6) 
fi^ = MG(l-cos0)/r^ 

where </> and are the zenith angle and the azimuthal angle in 
spherical coordinates centered at the Galactic center, respec- 
tively. /iG = 184.2/iG kpc-' is the magnetic moment of the 
Galactic dipole, which is normalized at 0.3/iG in the vicinity of 
the solar system. Theoretically, a dipole field c an be predicted 
by the so-ca lled AO mode in the dynamo theory jWidrow 120021; 
lBeck1l20 09'). Since the AO mode predicts the A-type parity of a 
spiral component, it is non-trivial that GMF models with S-type 
parity have a dipole component naturally. Note that there is no 
clear observational evidence of a dipole magnetic field. 

GMF has also a significant turbulent component, but its field 
strength, Btuib, has large uncertainty. Several observations have 
estimated as Btuib =(0.5 - 2) Z?reg (Beck 2001), where Bi-eg is the 
strength of the regular component of GMF. In this study, we as- 
sume Btuib = Bres and Bieg is the sum of the strength of the spiral 
component and the strength of the dipole component. The co- 
herent length of the turbulent component, Ic, has been thought 
to be typically < 100 pc (Beck 2001). We assume a turbu- 
lent magnetic field with the Kolmogorov power spectrum with 
Ic = 50 pc. Since this Ic is much smaller than the Larmor radius 
of HEPs the turbulent component does not significantly affect 
the pr opagation of HEPs with ~ 10^*' eV ( Tinvakov & Tkachev I 
l2005h . However, this component will be important for the prop- 
agation of heavy nuclei in the Galactic space due to their large 
charges. 

In this study, we adopt a tu rbulent IGMF model constructed 
in lYoshiguchi et al. I (l2003al) . We cover Galactic space with 
cubes of side Ic and assume that turbulent magnetic field in each 
cube is represented as a Gaussian random field with zero mean 
and a power-low spectrum 

(B(k)B*(k)) I -(J"" ^1?;/';^ ' ^ 2'^/'-' (7) 

where /cut is a numerical cutoff scale. Physically, one expects 
'cut ^ but we set Zcut = to save the CPU time. We adopt 
n-n=-\l /2>, which corresponds to the Kolmogorov power spec- 
trum. 

2.2. Source model of HEPs 

The PAO data implies that HECR sources are distributed 
to be compara ble with matter (or ga laxy) distribution in lo- 
cal Universe d Abraham et al. I l2007l [20 08; Ghiselli ni et al. I 
120081: iKashti & Waxman I l2008l 'Ta kamTet al. ii200 9). Follow- 
ing this result, we consider HEP source model in which sources 



are distributed following galaxy distribution actually observed 
Takami et al. (2006). This model is summarized as below. See 
OTakami et al. i i2006) for more details. 

This source model is based on the Infrared Astronomical 
Satellite Poi nt Source Redshift Su rvey (IRAS PSCz) catalog 
of galaxies dSaunders et al. 1 12000|) . The IRAS PSCz cata- 
log is appropriate to construct the model of HEP sources dis- 
tributed with being comparable to galaxy distribution because 
this galaxy catalog is a flux-limited, i.e., uniform, catalog of 
galaxies with very large sky coverage 84% on the whole 
sky). This catalog includes 2 selection effects: One is the ef- 
fect of undetected dark galaxies with the flux below the flux- 
limit (0.6 Jy), and the other is ~ 16% of unobserved area which 
mainly corresponds to the sky near the Galactic plane. We cor- 
rect these selec tion effects by us ing a luminosity function of the 
IRAS galaxies ( iTakeuchi et al. 1200 3). The former effect can be 
corrected by adding virtual galaxies with the flux less than the 
flux-limit in the vicinity of the original IRAS galaxies not to 
spoil the original structured galaxy distribution. For the latter 
effect, we simply assume that galaxies are distributed isotrop- 
ically in the unobserved area. These corrections allow us to 
construct (a model of) a complete galaxy catalog. We adopt this 
corrected galaxy catalog only within 200 Mpc to make HECR 
source distribution because the number of original galaxies out- 
side 200 Mpc is too small to make source distribution reflecting 
matter distribution in the Universe. It is assumed that galaxies 
are isotropic ally distributed outside 200 Mpc up to 1 Gpc. 

The subsets of this corrected galaxy catalog are regarded 
as HEP source distributions, i.e., we randomly select galaxies 
from this corrected catalog according to a given number density 
of HEP sources, n,, and then the selected galaxies are regarded 
as HEP sources. The distribution of HEP sources traces the dis- 
tribution of the original IRAS galaxies within 200 Mpc. HEP 
sources outside 200 Mpc are distributed isotropically by defi- 
nition, but these isotropic sources hardly contribute to the total 
flux of HEPs because more than 90% of HEPs above 6 x 10'^ 
eV detected at the Earth is generated within 200 Mpc. 

For the physical properties of HEP sources, we assume that 
all the sources emit HEPs persistently with the same power 
and with a power-law spectrum of oc E~^-^, where E is the 
energy of HEPs. This spectral index can well reproduce the 
observed spectrum of HECRs above 10''' eV. We consider 2 
number densities of HEP sources, «, = 10""^ and 10"^ Mpc"^. 
lO"'* Mpc""* can best reproduce anisotropy in the arrival distri- 
bution of HECRs d etected by the PAO (Takami & Sato 20^ 
ICuoco et al. 1 I2009I) . while 10"^ Mpc'^ is the number density 
to repr oduce the AGASA anisotropy well ("Takami et al. 1120061 : 
iBlasi & de Marco 2004; Kachelriess & Semikoz 2005^ 

When selecting HEP sources from the corrected galaxy cat- 
alog, we do not consider any specific condition, e.g., selecting 
only AGNs and so on. Thus, we generate 500 source distribu- 
tions for = 10"^ Mpc"^ and = 10"^ Mpc""*, and discuss the 
cross-correlation between HEPs and their sources statistically. 

2.3. Simulation of Arriving HEPs 

In order to simulate the arrival distribution of HEPs for a 
given so urce distribution, we adopt a calculation method devel- 
oped bv lTakami et al. I (12006). When a cosmic ray is injected 
into magnetized Universe from its source, whether this cosmic 
ray can reach the Earth or not cannot be known until its tra- 
jectory is actually calculated because of the existence of mag- 
netic field. Such a cosmic ray does not arrive at the Earth in 
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Fig. 1 . — Examples of trajectories of protons with tlie energy of lO"^ = 6.3 X lO" eV projected onto a plane perpendicular to the Galactic disk and including 
the Galactic center and the Earth. The position of the Earth is j: = 8.5 kpc. Only the spiral component of GMF is considered. The arrival directions of protons are 
(£,b) = (180.0°, 0.3°) {upper panel), (180.0°, -0.3°) {lower left panel), and (0.0° ,-44.8°) (lower right panel). The 4 trajectories in each panel correspond to the 
trajectories of protons in the BS-S model (red), BS-A model {green), AS-S model (blue), and AS-A model (magenta). 



many case. This requires the injection of many particles and 
the calculation of their trajectories to simulate the arrival dis- 
tribution of HECRs at the Earth. It is time-wasting to calculate 
the trajectories of HECRs which cannot reach the Earth. Our 
method solved this problem by considering only the trajecto- 
ries of HEPs reaching the Earth, and enables us to calculate the 
arrival distribution of protons within reasonable CPU time. 

In this method, the trajectories of protons are calculated 
by the backtracking method. The trajectory of an oppositely- 
charged proton ejected from the Earth can regarded as that of 
a proton coming from extragalactic space. Many oppositely- 
charged protons are ejected from the Earth isotropically, and 
their trajectories in the Galactic space are calculated. A decade 
of the energy of HEPs is divided into 10 bins in logarithmic 
scale. 125,000 particles are emitted from the Earth per energy 
bin in this study. We consider HEPs with the energy from 10'^ 
eV to 10^' eV at the Earth. Once the ejected particles escape 
from the Galactic space, they propagate straightforwardly in in- 
tergalactic space with energy-ga/n. The boundary between the 
Galactic and extragalactic space is set to be 40 kpc from the 
Galactic center The calculation of their trajectories is finished 
when the comoving distance of the particles exceeds 1,500 Mpc 
or the energy of the particles exceeds 10^^ eV. This energy 
corresponds to the maximum acceleration energy of protons at 
sources. As the energy-loss processes of HEPs (the energy-ga/n 
processes of oppositely-charged protons), in addition to Bethe- 
Heitler pair creation and photopion production by interactions 
with cosmic microwave background (CMB) photons, adiabatic 
energy-loss due to the cosmic expansion is also considered in 
intergalactic space. All the energy-loss can be neglected dur- 



ing propagation in the Galactic space because all the energy- 
loss lengths of protons are much larger than the scale of the 
Galactic space. We assume the standard ACDM cosmology 
with VL„ = 0.3, r^A = 0.7, //o = 71 km s"' Mpc"'. An analytical 
fitting form ula of the energy-lo s s leng th of Bethe-Heitler pair 
creation by IChodorowski et alTI d 1 9921) is used. For photopion 
production, we adopt inelasticity and mean-free path calculated 
by an event generator SOPHIA (iMucke et al. 11200 0*). 

When the trajectory of the ith oppositely-charged proton 
passes over sources, a weight below is attached into this tra- 
jectory for a positive arriving proton, 

P rFnn..-V ^ dN/dE,(E,') dE, 



dE 



where j labels sources on the trajectory. Zi.j and dj j are 
the redshift and comoving distance of the jth source, respec- 
tively. £'1,'' is the energy of the particle at the jth source. 
dN /dEg(Eg) oc Eg~^'^ is the injection spectrum of HEPs at the 
jth source. dEg/dE can be calculated by the propagation data. 
£■"' " is a factor which originates from the ejection spectrum 
of oppositely-charged protons at the Earth {dN / dlog^^E (x 
const.). This weight corresponds to a relative probability in or- 
der that the ith particle is a really arriving proton. This factor is 
attached to all the trajectories of oppositely-charged protons. In 
order to obtain protons arriving at the Earth, we randomly select 
the required number of trajectories according to the probability 
proportional to Pseiec(£^,0- The energy and arrival directions of 
the protons are the energy and the injection directions of the 
corresponding oppositely-charged protons. We assume that the 
arrival directions fluctuates following a 2 dimensional Gaussian 
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Fig. 2. — Arrival dii'ections of protons with the energy above lO''' eV emitted from 8 sources located at 50 Mpc from the Earth {black circles) in Galactic 
coordinates. The energy of each proton is represented in color The 4 different spiral components of GMF are taken into account. 



distribution with = 1 ° to take the finite accuracy of the arrival 
direction determination of HECR observatories into account. 

2.4. Statistical Methods 

We will discuss the correlation between the arrival di- 
rections of HEPs and their sources in the next section. 
A cross-correlation function for this purpose is defined as 
dTakami et al. l2009HBlakeetal. Il2006h . 



ES(9) - ES'iO) - E'S(0) + E'S'iO) 
E'S'(e) ' 



(9) 



where E and S represent HEP events and their sources which are 
considered to calculate the cross-correlation function, respec- 
tively. ES{0) is the number of pairs between E and S with the 
separation angles from 6* to 6* + divided by A^A^, for normal- 
ization, where and A'^, are the number of events and the num- 
ber of sources considered, respectively. A6' is set to be 1°. E' 
and S' represent randomly distributed events following the aper- 
ture of a HECR observatory and randomly distributed sources 
following the selection window of a source catalog, respec- 
tively. Here, S' represents just sources randomly distributed in 
the sky because the selection effects were fixed. ES'(O), E'S{9), 
and E'S'{9) are defined following the definition of ES{9). E' 
and S' allow to correct their non-uniform apertures. By defini- 
tion, Wes{9) > means positive correlation. 

The aperture of a ground array is simply described as 
dSommers 2001) 

uj(S) oc cos(ao) cos((5) sin(am) + a,„ sin(flo) sin((5) , (10) 

where a,,, is given by 





TT 



'(0 



if.e> 1 

if.e<-l 

otherwise 



and 



cos(8) - sin(flo) sin{5) 
cos(flo)cos(5) 



(11) 



(12) 



Here, oq and Q are the terrestrial latitude of a ground array 
and the zenith angle for an experimental cut. We simulate the 
apertures of 2 HECR observatories; the PAO and Telescope 
Array (TA). These 2 observatories are complementary because 
the PAO observes the southern sky while the TA observes the 
northern sky. Thus, we consider these 2 apertures as repre- 
sentatives of HECR observatories in the southern and north- 
ern sky. Parameters are aq = -35.2° and O = 60° for the PAO 
(lAbraham et al 1120081) . and ao = 39.3° and 6 = 45° for the TA 
dNonaka et an2009l) . 



3. RESULTS 

3.1. Deflections of protons by GMF 

Before considering the correlation between the arrival direc- 
tions of HEPs and their sources, we examine the trajectories of 
HEPs in OMR 

The dipole component of GMF is directed toward the north 
Galactic pole in the Galactic plane, and therefore protons com- 
ing from extragalactic space propagate clockwise, as viewed 
from the north Galactic pole. This fact indicates the arrival di- 
rections of HEPs are shifted from the positions of their sources 
to lower Galactic longitude. This effect is weaker in the di- 
rection of the Galactic anti-center because of the smaller field 
strength than in the direction of the Galactic center The strong 
dipole field near the Galactic center makes the complex trajec- 
tories of HEPs. In addition, since the dipole field above and 
below the Galactic center has the opposite direction to that in 
the Galactic plane away from the Galactic center, the trajec- 
tories of HEPs from the direction of ^ ^ 0°, where £ is the 
Galactic longi tude, can be defl ected anticloc kwise in the direc - 
tion of ^ ~ 0° dYoshiguchi et a l. 2004; Taka mi & Sato l2008bl) . 
The deflection angles of HEPs are not so large that the cor- 
relation disappears completely (< 10°) except the direction of 
the Galactic center if p rotons above 10'^ eV are considered 
dTakami & Sato Il200 8b^. 

The turbulent component of GMF makes HEPs propagate 
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Fig. 3. — Distribution of the deflection angles of protons with the energy of 10 eV by the spiral component of GMF as a function of the arrival directions of 
them in equatorial coordinates. Their deflection angles are shown in color. The difference between the AS models and BS GMF models is prominent especially in 
the northern hemisphere. 



diffusively. The typical deflection angle of a HEP is 

0(e„DGal)^l°(^^) (^), 

(13) 

where ep and Dcai are the energy of protons and the typical 
size of the Galaxy, respectively. This estimation is expected 
to be larger than the deflection angles of protons by the turbu- 
lent field in our simulation, since the turbulent field becomes 
weaker outside the Galactic disk. Thus, the turbulent compo- 
nent of GMF subdominantly contribu tes to the total deflection 
of HEPs dTinvakov & Tkachev II200I and the global deflection 
pattern of HEPs produced by t he coherent comp onents of GMF 
is not largely disturbed ( Yoshig uchi et al. Il2004l) . 

The spiral component of GMF dominantly affects the deflec- 
tions of HEPs. This is because the solar system is located in 
the Galactic disk, where the spiral field is stronger than in the 
Galactic halo. The spiral field is stronger (^ l.SfiG) than the 
dipole field (~ 0.3/iG) in the vicinity of the solar system. Sev- 
eral examples of the trajectories of HEPs with 10'^ '^ eV, which 
are useful to understand the features of the trajectories in the 
spiral field, are shown in Fig. [1] Only the spiral component is 
taken into account as GMF. The upper panel shows the trajec- 
tories of HEPs in the 4 different spiral field models in x-z plane. 
Note that ;ic-axis is a line penetrating from the Galactic center to 
the solar system, and z-axis is a line penetrating from the Galac- 
tic center to the north Galactic pole. The solar system is located 
at X = 8.5 kpc. The arrival directions of these HEPs at the Earth 
are the same, (i,b) (180.0°, 0.3°). Practically, we calculate 
the trajectories of oppositely-charged protons injected from the 
Earth to i£,b) ^ (180.0°, 0.3°) and regard these trajectories as 
those of arriving HEPs. The trajectory in the BS-S model and 
that in the BS-A model are duplicated because these two models 
has the same spiral field above the Galactic plane. The trajec- 
tory in the AS-S model and that in the AS-A model are in the 
same situation. The trajectories in the BS models are separated 
from those in the AS models outside 15 kpc because of the 



existence of the field reversal of the spiral field in the BS mod- 
els. Once a proton passes over the radius of a field reversal, the 
deflection direction of the proton is changed into the opposite 
direction to the direction before. This effect does not occur in 
the AS models, since the AS models do not have the field rever- 
sals of GMF. This effect also allows the total deflection angles 
of HEPs to be smaller in the BS models than in the AS models. 
After the oppositely-charged protons propagate away from the 
Galactic plane, their deflections are smaller because the spiral 
component of GMF becomes smaller exponentially. Thus, the 
spiral GMF in the vicinity of the solar system mainly affects the 
total deflections of HEPs. 

The effect of the parity of the spiral component of GMF can 
be observed in the other panels. The lower left panel is the 
same as the upper panel, but the arrival directions of HEPs are 
ie,b) = (180.0°, -0.3°). Due to the difference of the parity, the 
spiral field of the BS-A (AS-A) model has the opposite direc- 
tion to that of the BS-S (AS-S) model below the Galactic plane. 
In this panel, the trajectories of oppositely-charged protons are 
deflected toward the upside of the Galactic plane near the so- 
lar system in the cases of GMF models with the S-type parity 
(BS-S and AS-S), while the trajectories are deflected toward 
the downside of the Galactic plane for the GMF models with 
the A-type parity (BS-A and AS-A). Compared to the upper 
panel, the BS-A and AS-A models show the symmetric deflec- 
tion pattern. The lower right panel is another example to show 
the effect of the parity. HEP trajectories with the same arrival 
direction, (e,b) = (0.0°, -44. 8°), are separated by the effects of 
the parity and field reversals of GMF. 

The deflection directions of HEPs depend on the positions 
of their sources. Fig. |2] demonstrates the arrival directions of 
protons emitted from 8 extragalactic sources at 50 Mpc from 
the Earth. The positions of the sources are shown by black cir- 
cles. Only the spiral component of GMF is taken into account. 
The positions of the sources are artificially set and protons with 
the energy down to 10'^ eV are considered to see the directions 



Does GMF Disturb the Correlation of HECRs with their Sources? 



7 



15 



10 



rij^lO"' Mpc"-', Np = 27 
E > 60 EeV, d < 75 Mpc 



No GMF 
BS-S model 
AS-A model 



- 



5 10 15 

Separation Angle: 6 [deg] 



20 



120 
100 
80 
60 
40 
20 

-20 



n5=10"^ Mpc"^, Np = 27 
E > 60 EeV, d < 75 Mpc 



No GMF 
BS-S model 
AS-A model 



5 10 15 

Separation Angle: 6 [deg] 



20 



Fig. 4. — Cross-correlation functions between the arrival directions of 27 simulated HEPs with the energy above 6 X lO" eV and the positions of HEP sources 
with ris ~ 10^ (left) and Kr'Mpc"^ (right) within 75 Mpc. The aperture of the PAO is taken into account. Only the spiral component is considered as GMF; no 
GMF (green), the BS-S model (red), and the AS-A model (blue). 



and strengths of the deflections of HEPs. The energy of protons 
are represented by color in the figure. In order to make these 
panels, we calculate the trajectories of oppositely-charged pro- 
tons injected from the Earth, select trajectories which pass the 
sources, and plot the injection directions of the particles. 

The arrival directions of HEPs are basically shifted from 
the positions of their sources to the directions of low 
or high latitude dStarievI 1 19971: lAlvarez-Muniz et al. 1 120021: 
iTakami &Satoll2008'bir and are arranged reflecting their ener- 
gies beca use the spiral field is dire cted parallel to the Galac- 
tic plane dYoshiguchi et al. Il200 3b'). The deflection directions 
depend on the structure of the spiral field. Protons emitted 
from a source with (£,b) ^ (120°, 60°) are deflected to lower 
latitude. We can see that the deflection angles of these pro- 
tons are smaller in the BS-S model and BS-A model than in 
the AS-S model and AS-A model, as discussed above. Since 
the AS-A model and BS-A model have symmetry above and 
below the Galactic plane, the arrival directions of HEPs from 
the source and those from a source with i£,b) ^ (120°, -60°) 
are symmetric. This symmetric feature can be seen for HEPs 
emitted from the other sources. Comparing protons from a 
source with i£,b) ~ (120°, 60°) with those from a source with 
(£, b) ~ (-60°, 60°), for example, we can see the dependence of 
the deflection angles of protons on source position. The latter 
protons are less deflected even for ^ lO''''^ eV than the for- 
mer protons. We can also see that protons with the energy of 
^ 10'^ eV are trapped near the Galactic center and are largely 
scattered. 

The arrangement of the deflection directions of protons can 
be a hint to un derstand the structure of GMF in the direc - 
tion of sources dYoshiguchi et al. Il2003bt iGolup et al. Il2009h . 
In fact, protons with the energy lower than ~ 4 x 10'^ eV 
are difficult to be a probe of GMF structure, since cosmolog- 
ically distant sources can significantly contribute to the total 
flux of HECRs. The background makes the identification of 
the sources of such low energy protons difficult. On the other 
hand, HEPs above - 6 x lO''' eV is a good tool for this pur- 
pose because they can reach the Earth only from nearby sources 
due to Greisen-Zatsepin-Ku z'min (GZK) effect dGreisenll966l: 
IZatsepin & Kuz'min |[T966l) . 

Finally, we show how the deflection angles of protons are dis- 
tributed in the sky as a function of their arrival directions. Fig. 
[3]shows the distribution of deflection angles of protons with the 
energy of 10'^^ eV in equatorial coordinates. Only the spiral 



component of GMF is taken into account. In all the models, 
the trajectories of HEPs are highly deflected near the Galactic 
center. We can also see that the BS models predict the smaller 
deflection angles of HEPs than the AS models, as discussed in 
Figs. [T]and|2] The difference between the AS models and BS 
models is noticeable especially in the northern sky, where the 
Galactic center is not contained. A similar plot taking the d ipole 
field into account can be seen in iTakami & Satol d2008bl) . The 
effect of the dipole field appears essentially only in the direction 
of the Galactic center 

3.2. Cross-correlation between HEPs and their sources 

First of all, we discuss the effect of the spiral component of 
GMF on the cross-correlation function between the arrival di- 
rections of HEPs and the positions of their sources. The effects 
of the dipole and turbulent components of GMF will be consid- 
ered afterwards. Fig. |4] shows the cross-correlation functions 
between the arrival directions of simulated HEPs with the en- 
ergy above 6 x 10'^ eV and the positions of their sources with 
n.5 ^ 10""* {left) and 10"^Mpc"^ {right) within 75 Mpc in simu- 
lation. The number of HEPs is set to be 27 events, which is the 
same as that of the published PAO events. The aperture of the 
PAO is taken into account. These cross-correlation functions 
are calculated under 3 different GMF models, i.e., no GMF 
{red), the BS-S model {green), and the AS-A model {blue). The 
points represent the averages of the cross-correlation functions 
over 500 source distributions and the error bars correspond to 
1(7 errors, which means that 68% of the values of the cross- 
correlation functions over the 500 realizations is contained in 
the range of the error bars. 

The averages and errors are estimated as follows: we sim- 
ulate one arrival distribution of a required number of HEPs 
(27 in this case) for a given source distribution, and then cal- 
culate the cross-correlation function between simulated events 
and sources within 75 Mpc in this simulation. We repeatedly 
calculate the cross-correlation functions in the same way for 
500 different source distributions with the same n,. Then, we 
calculate the averages of these 500 cross-correlation functions, 
and the values of the cross-correlation functions at 16% from 
the minimum and at 84% from the maximum of these func- 
tions. The latter 2 values are the lower limit and the upper limit 
of the error bars, respectively. Therefore, the error bars include 
not only errors due to the finite number of simulated events but 
also errors originating from the uncertainty of HEP source po- 
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Fig. 5. — Cross-correlation functions of 27 simulated HEPs with the energy above 6 X lO" eV and their sources with «j = 10"* Mpc-^ within 75 Mpc. The BS-S 
model (red), the BS-A model (green), the AS-S model (blue), and the AS-A model (magenta) are considered as a spiral component of GMF. The apertures of the 
PAO (left) and TA (right) are taken into account. 



sitions, i.e., the galaxy sampling in our HECR source model 
(see ^2.2b . Note that the latter error cannot be reduced unless 
additional conditions are required in the source model. All the 
cross-correlation functions except for Fig. |9]are plotted in the 
same way as above. 

The cross-correlation functions have a sharp peak at the 
smallest angular bin for no GMF in the both panels of Fig. |4] 
because the propagation trajectories of HEPs are straight. The 
finite accuracy of HECR observatories to determine the arrival 
directions of HEPs 1°) allows the positive values of the 
cross-correlation functions at the second and third bins. The 
value of the cross-correlation function at the peak is larger for 
ris = 10"^ Mpc"^ than that for n, = 10"^ Mpc"-' because a smaller 
number of events is enough to trace the positions of sources by 
smaller number density. When the BS-S model or the AS-A 
model of GMF is considered, the positive signals of the cross- 
correlation functions for no GMF are suppressed. The cross- 
correlation functions for both models are consistent with no 
correlation (Wes(6') = 0) within la error in both panels except 
at several angular bins for the BS-S model. 

The left panel of Fig. |5]is the same as the left panel of Fig. 
m but the cross-correlation functions for the BS-A model and 
the AS-S model are added and the result for no GMF is re- 
moved. The cross-correlation functions for the AS models are 
consistent with no correlation, while those for the BS models 
have a signal of positive correlation at 3°-4°, but significance 
is low. This results from the fact that the BS models predict 
smaller deflection of HEPs than the AS models. It is suggestive 
that the angular s eparation scale of the PAO correlation is 3.1° 
dAbraham et al. | [2Q07, 2008). However, note that the analysis 
of the PAO was based on cumulative cross-correlation, which is 
different from differential cross-correlation used in this study. 

The right panel of Fig. |5]is the same as the left panel, but the 
aperture of the TA is taken into account to observe the differ- 
ence between a cross-correlation signal in the southern sky and 
that in the northern sky. The number of simulated events is the 
same as for the PAO, 27. Note that it spends more time than 
the PAO that the TA accumulate 27 events above ^ 6 x 10'^ eV 
because the aperture of the TA is smaller than that of the PAO. 
In the northern sky, the AS models predict the cross-correlation 
functions consistent with no correlation within 1 a error in all 
the angular scale. Since the lower limit of the error bars is 
lower than that for the left panel, detecting the correlation in 
the northern sky is more difficult than in the southern sky if ax- 



isymmetric GMF is realized. Comparing both panels of Fig. |5] 
we notice that the cross-correlation functions predicted by the 
BS-S (AS-S) model is quite similar to those predicted by the 
BS-A (AS-A). The parity of the spiral component of GMF does 
not largely affect this statistical quantity, although the parity 
produce the difference of the deflection directions of HEPs, as 
seen in Fig. |2l The difference between the AS models and the 
BS models is more clear in the northern sky, as we will show 
afterwards. 

The error bars of the cross-correlation function originate 
from not only the finite number of events but also the uncer- 
tainty of positions of HECR sources (i.e., the galaxy sampling 
in our source model). The former errors can be decreased by 
increasing the number of detected events, while the latter errors 
cannot be reduced unless additional constraints are given in our 
model. For an ideal case where the number of detected events is 
infinite, we can consider only the latter errors. In the practical 
point of view, it is meaningful to estimate the number of events 
required for saturation of the errors because we can discuss the 
correlation of HEPs and their sources based on almost only the 
uncertainty originating from the galaxy sampling after accumu- 
lating such number of events. Fig. |6]shows the cross-correlation 
functions between HEPs with the energy above 6 x 10'^ eV and 
their sources with = lO"** Mpc"-' {left) and = 10"^ Mpc"-' 
(right) within 75 Mpc. The HEP events are simulated by as- 
suming the BS-S model and the PAO aperture. The numbers 
of protons considered are 27 (red), 100 (green), 200 (blue), 500 
(magenta), and 1000 (light blue). When the number of events is 
small, errors due to the finite number of events are dominated 
in the total errors. As increasing the number of events, these 
errors decrease and errors due to the sampling of galaxies be- 
come dominant. Since the lengths of error bars for 500 events 
are not largely different from those for 1000 events, the error 
bars are well saturated at ^ 500 events for both n,. This fact is 
unchanged for the other GMF models and in the northern sky. 
The PAO can accumulate ^ 200 events abov e ^ 6 x 10'^ eV 
for 10 years ( 27 events per 1.5 years exposure dAbraham et al. I 
2008)). The error bars for 200 events are also sufficiently small. 
Thus, we consider 200 protons for discussions below. 

Fig. |7] shows the cross-correlation functions between 200 
simulated HEPs with the energy above 6 x 10'^ eV and their 
sources with n, = 10"'* Mpc"^ within 75 Mpc for the BS-S 
model (red), BS-A model (green), AS-S model (blue), and AS- 
A model (magenta). The apertures of the PAO (left) and TA 



Does GMF Disturb the Correlation of HECRs with their Sources? 



9 




5 10 15 

Separation Angle: 9 [deg] 



20 



25 
20 
15 
10 
5 




= 10"^^ Mpc"^ 
E > 60 EeV, d < 75 Mpc 



-4. 



27 events 
1 00 events 
200 events 
500 events 
1 000 events 



5 10 15 

Separation Angle: 9 [deg] 



20 



Fig. 6. — Cross-correlation functions between HEPs simulated in the BS-S model with the energy above 6 X lO" eV and their sources with ris = lO"* Mpc"^ 
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Fig. 9. — Cross-correlation functions between 200 protons simulated from a given source distribution with ji, ~ 10 Mpc and their sources within 75 Mpc. 
The energy of the protons is above 6 X lO" eV. The apertures of the PAO (left) and TA (right) are taken into account. Only the spiral component of GMF is taken 
into account. 



(right) are taken into account. Compared with Fig. |5] the error 
bars become smaller, and the clear signals of positive correla- 
tion appear in some angular scale. In addition, we can see that 
the cross-correlation functions for the BS-S (AS-S) model be- 
have similarly to those for the BS-A (AS-A) model. In the left 
panel, the cross-correlation functions for the BS models are not 
consistent with zero (no correlation) up to 7° at Icr error, and 
therefore the positive correlation appears in this angular scale. 
On the other hand, the signal of positive correlation for the AS 
models is in the angular scale from 2° to 8°. The lower limit 
of the error bars at small angular scale is lower than that for 
the BS models. This reflects the fact that the AS models pre- 
dict the larger deflection of HEPs than the BS models. The 
right panel shows qualitatively the same result as the left panel. 
The BS models predict the positive correlation up to 6°, while 
the cross-correlation functions for the AS models are incon- 
sistent with zero at 3° -8°. The difference between the cross- 
correlation functions for the BS models and those for the AS 
models is larger in the northern sky than in the southern sky. 
Fig. |8]shows the same figure as for Fig. |7] but «, = 10"^ Mpc"^. 
The qualitative features of the cross-correlation functions are 
the same as those in Fig. |7] The values of the cross-correlation 
functions are larger than those for = 10"^ Mpc"^ because of 
the smaller n,, but the error bars also become larger 

The difference between the BS models and the AS models 
can be clearly seen in the cross-correlation functions calcu- 
lated from a given source distribution. In this case, the error 
bars of the cross-correlation functions are determined by only 
errors due to the finite number of events, i.e., the error bars 
can be minimized by accumulating events. Fig. |9] shows the 
cross-correlation functions based on one specific source distri- 
bution with ris = 10"^ Mpc"^ The apertures of the PAO (left) 
and TA (right) are considered. The arrival distribution of 200 
HEPs with the energy above 6 x 10'^ eV is realized 500 times, 
and then the averages and the 68% error bars of the 500 cross- 
correlation functions are calculated. For this example, the dif- 
ference between the cross-correlation functions for the BS mod- 
els and those for the AS models is clear in the northern sky. The 
BS models predict the strongly positive correlation within 5°, 
while the cross-correlation functions are consistent with zero at 
small angular scale for the AS models. On the other hand, the 
positive correlation is predicted for all the spiral field models 
in the southern sky. The angular scale at which the positive 



correlation appears is larger for the AS models than for the BS 
models. 

Next, we investigate the effect of the dipole component and 
the turbulent component of GMF on the cross-correlation func- 
tion. Fig. [TO]shows the cross-correlation functions between the 
arrival directions of 200 HEPs with the energy above 6 x 10^^ 
eV and the positions of their sources with «s = 10"^ Mpc"-' 
within 75 Mpc. The apertures of the PAO and the TA are taken 
into account in the left panel and the right panel, respectively. 
We consider the BS-S model in the upper panel and the AS-A 
model in the lower panel as a spiral component of GMF (red). 
In addition to the spiral component, the dipole component of 
GMF (green), the turbulent component of GMF (blue), and both 
the dipole and turbulent components (magenta) are considered 
for comparison. In all the panel, the dipole field reduces the av- 
erages and the positions of the error bars of the cross-correlation 
functions, and therefore the positive signals of the correlation 
decreases at small angular scale. On the other hand, the tur- 
bulent component of GMF hardly affects the cross-correlation 
functions. 

Finally, we check the dependence of the cross-correlation 
function on the maximum distance of HEP sources to calculate 
the cross-correlation functions. So far, the maximum distance 
has been fixed at 75 Mpc, which is the same dis tance as the 
analysis of the PAO (Abraha m et 111120071 12008|) ( z = 0.018 
in the concordance cosmology ). However, this value does 
not have inevitability, although it is consistent with the GZK 
scenario, since this value was derived to optimize the correla- 
tion signal. Thus, we check the cross-correlation functions for 
d < 100 Mpc and d < 50 Mpc. 

Fig. [TT] shows the cross-correlation functions of 200 HEPs 
with the energy above 6 x lO'^ eV and their sources with 
«.s ~ 10""* Mpc""* (upper) and n, ~ 10"^ Mpc"-' (lower) within 
100 Mpc. The apertures of the PAO (left) and TA (right) ai-e 
taken into account. Compared to Figs. |7]and[8] the qualitative 
features of the cross-correlation signals are the same. How- 
ever, the values of the averages and the positions of the error 
bars become smaller than those in Figs. |7]and[8] The lengths 
of the error bars also become smaller These are because the 
number of sources considered to calculate the cross-correlation 
functions increases due to increasing the maximum distance of 
the sources for calculating the cross-correlation functions. The 
large number of sources decreases the fluctuation of galaxy 
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sampHng. Fig. fT2l shows the same as Fig. (TT] but the maxi- 
mum distance of HEP sources considered to calculate the cross- 
correlation functions is 50 Mpc. These cross-correlation func- 
tions also have the same features qualitatively as those in Fig. 
nn Contrary to Fig. [TT] the values of the cross-correlation 
functions are larger than those in Figs. |7]and|8]because of the 
smaller number of sources considered to calculate the cross- 
correlation functions. The small number of sources also pro- 
duces the large fluctuations of the cross-correlation functions. 
Thus, the qualitative features of the positive correlation are not 
dependent of if t/ < 100 Mpc. 

3.3. Angular Distance Sensitive to Positive Correlation 



In ^3.21 we investigated the appearance of the positive cor- 
relation between HEPs and their sources based on a general 
source model in which HEP sources are distributed with being 
comparable to galaxy distribution. The consideration of 500 
source distributions allows us to conclude that the positive cor- 
relation can appear at angular separation scale less than 10° 
even considering the effect of GMF, but this angular scale de- 
pends on the modelling of GMF. On the other hand, the posi- 
tive correlation between HECRs and their source candidates is 
generally tested by comparing the experimental data of HECR 
events to ran dom distribution assuming one source candidates 
distribution Abraham et al. ' '2007", '2008'; 'Kashti & WaxmanI 
2008|;lGeorg e et al. 2008; GhiselHni et al. 2008; Takami et al. 
200% . In the analysis of the PAO (A braham et al. 120071,120081) . 
the energy threshold of HECRs, the maximum distance of 
source candidates to investigate the correlation, and the angular 
separation scale between the arrival directions of HECRs and 
the positions of source candidates were regarded as parameters. 



and the parameters were constrained so that the significance of 
excess events around source candidates over random event dis- 
tribution was maximized. (In fact, the PAO constrained these 
parameters by the first half of detected events, and then tested 
the correlation under the same parameter est by using the sec- 
ond half of detected events.) In this section, we estimate the 
angular scale at which the significance of the positive correla- 
tion between HECRs and their sources over random event dis- 
tribution is maximized. The energy threshold of HECRs and 
the maximum distance of sources to investigate the correlation 
are fixed to 6 x 10'^ eV and 75 Mpc, respectively. 

For random event distribution, the number of events which 
falls into circles with the angular radii 9, j, follows binomial 
distribution. 



Bi{N,pUOyj) = 



{Pr.n(0)y{^-{pUO)}f~\ (14) 



where is the total number of events. Pmn{9) is equivalent with 
the fraction of the sky (weighted by the exposure) defined the 
regions at angular separation les s than the angular sca l e from 
the selected source candidates (lAbraham et al. I l2007l 1200^ . 
Even when the effect of GMF is taken into account, we con- 
firmed that the number of events within 9 from the source can- 
didates follows binomial distribution. The effect of GMF is 
included in a parameter /7sim(^'), but Psim{9) is no longer the 
fraction of the sky regions. 

The most sensitive angular scale to the positive correlation 
can be found as follows. A source distribution is considered. 
First of all, we simulate (A^ =) 200 HEP events and count the 
number of events within 9 from the positions of sources, which 
is written as Msim- We repeat this step many times (10000 times 
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Fig. 1 1. — Cross-correlation functions of 200 simulated HEPs with the energy above 6 X lO" eV and the positions of their sources with ;ii ~ 10~* Mpc"^ 
(upper) and «j ~ 10"^ Mpc"' (lower) within 100 Mpc. The BS-S model (red), the BS-A model (green), the AS-S model (blue), and the AS-A model (magenta) are 
considered as a spiral component of GMF. The apertures of the PAO (left) and TA (right) are taken into account. 



in this study) and make the histogram of «sim- Then, we eval- 
uate Psim(^) by fitting this histogram by binomial distribution, 
Bi(2QQ,psim(Sy,j)- We can also estimate p^nid) from random 
distribution. The number of events inside circles centered at 
the source positions with radii less than 9 is n^^^- The proba- 
bility that the positive correlation is realized corresponds to the 
probability that h = nsim - "ran is positive. The probability dis- 
tribution following h is calculated by the convolution of the two 
binomial distributions. 



{Ps^(0)Y{\- 



E 



h + m 



Psim' 

N 
m 



(0)pue) 



m < N 
Since 



- n for 
this is 



(15) 



> 
func- 



at which " 

Etf-NP(N:P.n,(0),Pr 



we can estimate a 
is minimized, 6'be,st- 



(1-Psim(e))(l-Pran(^)) 

where m is summed over < 
and -n < m < N for n < 0. 
tion of 6 through Psim(^) and Prdn{0), 
.NP(N,P.iUO),pUdy,ii) 
\iy ,i^simi()),Pmn(S);n) gives the probability that the 
number of HEP events correlating their sources does not exceed 
the number of random events correlating the source positions. 
The values of p^imiS) and PrwiS) depend on source distribution 
because the apertures of experiments are not uniform and the 
effect of GMF depends on the arrival directions of HEPs. Thus, 
we calculate 6'best for 500 source distributions for each n,. 

The left panels of Fig. [13] shows the distributions of 6'best 
for «, = 10"^ Mpc"^ (upper) and 10"^ Mpc"^ (lower). The BS- 
S and AS-A models are considered as a spiral component of 
GMF. The dipole component and turbulent component are also 
taken into account. The apertures of the PAO (South) and the 



TA (North) are taken into account. The peak of 6'best distribu- 
tions strongly depend on the spiral shape of GMF. In the case of 
«.s = 10"^ Mpc"-', the peak of 6'best distribution is ^ 5° for the BS- 
S model, while that is 6° - 8° for the AS-A model. The cases of 
n.5 = 10"^ Mpc"^ predict larger 6'best at the peak; ^ 6° for the BS- 
S model and 8° - 10° for the AS-A model. The peak positions 
are smaller for n, = 10""* Mpc"^ than that for = 10"^ Mpc"^. If 
the number density of sources is large, source distribution pro- 
jected onto the sky is dense. In this case, the source closest to 
the arrival direction of a HEP is sometimes not a real source. 
Thus, the separation angles between HECRs and sources tend 
to be underestimated for larger n,. The difference between 0best 
distribution in the northern sky and that in the southern sky can 
be seen for the AS-A model. The peak position in the southern 
sky is ^ 2° smaller than that in the northern sky, although the 2 
distributions overlap. 

The right panels of Fig. [13] shows the distribution 
of the probability that n < by chance at 6'best, i-C-, 

I]fc-/v^(200,Psim(6'be.st),A-an(^'best);«) Parameters and GMF 
models used are the same as those in the corresponding left 
panels. In both n,, the BS-S model shows smaller chance prob- 
ability than the AS-A model. Almost all the source distribu- 
tions predict the chance probability less than 10"^, which cor- 
responds to the significance of the positive correlation more 
than ^ 3(7 in the context of Gaussian distribution, for the BS-S 
model. On the other hand, the AS-A model can predict larger 
chance probability because of the larger deflection of HEPs. 
Less than 5% of source distributions predicts chance proba- 
bility larger than 10"^ for n., = 10"^ Mpc""* . The situation is 
worse for n , = 10""* Mpc"^ because of the larger number density. 
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Fig. 12. — Same as Fig. 1111 but for HEP sources within 50 Mpc. 
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Fig. 13. — (left): Frequency distributions of the separation angle between 200 simulated HEPs with the energy above 6 X lO" eV and their sources with n., = 10^ 
IVIpc"^ (upper) and «, = 10"^ Mpc"^ (lower) within 75 IVIpc which gives maximal significance against random event distribution, 9bcsi> over 500 source distributions. 
The apertures of the PAO (South) and the TA (North) are taken into account. The BS-S model and AS-A model are considered as a spiral component of GMF. 
The dipole component and turbulent component are also taken into account, (right): Frequency distribution of the probability that the number of simulated events 
inside circles with the radii of 9\^^t centered at the position of their sources does not exceed to that for randomly distributed events. In other words, these are the 
distributions of the probability that the positive correlation does not occur All the parameters used are the same as those in the corresponding left panels. 
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Fig. 14. — (left): Frequency distributions of tlie separation angle between 200 simulated HHPs with the energy above 6 X 10 eV and their sources with ris = 10 
Mpc"^ within 75 Mpc which gives maximal significance against random event distribution, Sbcsti "ver 500 source distributions. The aperture of the PAO is taken 
into account. The BS-S model (upper) and AS-A model (lower) are considered as a spiral component of GMF (red). In addition to the spiral component, the dipole 
component of GMF (green), the turbulent component of GMF (blue), and both the dipole and turbulent components of GMF (magenta) are taken into account. 
(right): Frequency distribution of the probability that the number of simulated events inside circles with the radii of S(,est centered at the position of their sources 
does not exceed to that for randomly distributed events. In other words, these are the distribution of the probability that the positive correlation does not occur. All 
the parameters used are the same as those in the corresponding left panels. 



Does GMF Disturb the Correlation of HECRs with their Sources? 



15 



About 50% of source distributions predicts chance probability 
larger than 10"^. In other words, about 50% of source distri- 
butions does not produce the positive correlation at the signifi- 
cance more than 3ct. 

The dependence of these values on the dipole component and 
turbulent component of GMF is seen in Fig. [141 «j = 10""* 
Mpc"^ is considered in this figure. The left panels of Fig. [14] 
show 0best distributions for the BS-S model (upper) and the AS- 
A model (lower). In each panel, we consider the 4 configura- 
tions of GMF, only a spiral component, spiral component plus 
dipole component, spiral component plus turbulent component, 
and all the 3 components. The turbulent component hardly af- 
fects 6'be,st distribution contrary to the dipole component. The 
dipole component moves the ^^best distribution to the right side 
by - 0.5° for the BS-S model and - 1° for the AS-A model. 
Thus, additional regular components of GMF to the spiral field 
can significantly affect the angular scale of positive correlation. 
This effect of the dipole field is also observed in the right pan- 
els of Fig. [14] The dipole component of GMF increases chance 
probability. If the dipole field does not exist, only 5% of source 
distributions produce chance probability larger than 10"^. 

One often uses a penalty factor to estimate real sig- 
nifica nce when the experimental data of HECRs is ana- 
lyzed dTinvakov & Tkachev 1 1200 IL 120041: [Finlev & Westerhoff I 
|2004|) . The small number of HECRs detected at present indi- 
cates low significance of the correlation even if the correlation 
is true. Therefore, the significance must be accurately estimated 
to confirm the correlation. In order to estimate the significance, 
one scans over parameters (e.g., the maximum distance (red- 
shift) of sources, the energy threshold of HECRs, angular sep- 
aration, and so on) and finds the parameter set which gives the 
largest significance of the positive correlation. However, this 
significance is not real significance because this is optimized to 
experimental data, which is only one realization of HECR gen- 
eration. Thus, one estimates a penalty factor originating from 
the optimization of the parameter sets and corrects the signifi- 
cance. The need of the penalty factor originates from the fact 
that experimental data is only one realization. However, we can 
make many realizations of HEP events by simulation in this 
study and can estimate a real significance, as performed above. 
Thus, we do not need to consider a penalty factor in this study. 

4. DISCUSSION & CONCLUSION 

In this study, we considered the effect of GMF on the ex- 
pected positive correlation between the arrival distribution of 
HEPs with the energy above 6 x 10'^ eV and the positions of 
their sources by using simulations. We found that the positive 
correlation can be detected at angular separation scale less than 
10° even considering plausible GMF models when 200 HEPs 
above 6 x 10^^ eV are detected. In addition, the dependence 
of the positive correlation on GMF models was examined. The 
trajectories of HEPs are deflected more largely in the AS mod- 
els than in the BS models because of the lack of field rever- 
sals. This fact shows that the signal of the positive correlation 
is weaker for the AS models than for the BS models. The differ- 
ence between the AS models and the BS models is prominent 
especially in the northern sky. The parity of the spiral field 
does not affect the correlation signals. The dipole field, which 
is a model of the vertical component of GMF, can subdom- 
inantly contribute to the total deflections of HEPs, while the 
turbulent component of GMF hardly affects the deflections of 
HEP trajectories. Furthermore, we estimated the angular scale 



at which the significance of the positive correlation over ran- 
dom event distribution is maximized. The most sensitive angu- 
lar scale to the positive correlation depends on the spiral shape 
of GMF; 6'be,st ^ 5° for the BS models and - 7° for the AS mod- 
els. The significance of these positive correlation well exceed 
^ 3(7 in almost all the source distributions with n., = 10""* and 
10"^ Mpc""*, but 50% of source distributions with = 10"^ 
Mpc"^ predicts the positive correlation at the significance less 
than 3(7. This problem can be avoided by neglecting the dipole 
component. The dipole component contributes to 6'best by 0.5° 
- 1°. 

We assumed a dipole field as a vertical component of GMF, 
but the existence of a dipole magnetic field has not been not 
confirmed. There is large uncertainty on a vertical compo- 
nent of GMF at present. Galactic wind can be also a source 
of the vertical component of GMF, can predict vertical mag- 
netic field different from the dipole shape (Brandenburg et al. I 
Il993b . Recent soft X-ray observations h ave indicated the ex - 
istence of Galactic wind in the Galaxy dEverett et al. Il20()8h . 
Furthermore, high-sensitivity observations of several edge-on 
galaxies with galactic wind have revealed so-called X-shaped 
magnetic field in their halo (Krause 2007; Beck 2009). Al- 
though it is uncertain that the Galaxy has such non-dipole ver- 
tical component of GMF, such a component could contribute to 
the deflections of HEPs if exists. 

In addition to GMF, IGMF can also contribute to the total 
deflections of the trajectories of HEPs. However, since our 
current understanding on IGMF is poor, the deflections highly 
depend on IGMF modelling. A numerical simulation on cos- 
mological structure formation in local Universe showed a mag- 
netic struct ure is highly concentrated following matter distribu- 
tion ( Dolag et al. 1120051) . Since cross-sectional area occupied 
by strong IGMF is small, HEPs are deflected less than 1° ex- 
cept in the directions of clusters of galaxies. A simple mod- 
elling of IGMF by Takami et al. (2006) also predicted small 
deflection angles of HEPs and showed that HEP sources are 
unve iled by the arrival dire ctions of HEPs within a few de- 
gree dTakami & Satoll2008al) . These 2 IGMF models are based 
on realistic matter distribution of local Universe constructed 
by the IRAS catalog of galaxies ( SjiUJiderseial- 2000). On 
the other hand, hydro dynamical simula tion s on cosmologica l 
structure formation bv'S igl et alTI (|2004 and lDas et aTl (|2()08|) 
predicted relatively strong IGMF in filamentary structures up 
to lOnG. The large cross-sectional area occupied by filamen- 
tary structures with ^10 nG fields makes HEPs be deflected 
significantly. Note that these IGMF models have uncertainty 
on observer's positions because their resultant matter distribu- 
tions do not correspo nd to local Universe actually observed. 
The IGMF model of ISigl et al.1 ( |2004|) deflects the trajecto- 
ries of HEPs by ~ 20° on average even for 10^" eV. Thus, the 
large amount of information on the positions of their sources is 
lost. However, the correlation between HEPs and matter dis- 
tribution is, not alw ays l ost even under strong IGMF models. 
iKotera & Lemoine I (12008) pointed out the possibility of spuri- 
ous correlation based on another IGMF modelling. In this sce- 
nario, HECRs are scattered at strongly magnetized regions far 
from the positions of their sources and then reach the Earth. In 
this case, we can observe the correlation between H EPs and last 
scattering regions, which trace matter distribution. iRvu et al~l 
( |2009|) al so showed similar possibility under the IGMF mod- 
elling of iDas et al. ] (I2008h . As discussed above, the uncer- 
tainty of IGMF modellings does not allow us to conclude that 
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the correlation between the arrival directions of HEPs and their 
sources can be observed. The future projects of all-sky survey 
of Faraday rotation measurement, like Square Kilometer Array 
will give us useful information on IGMF in local Universe. 

Although we assumed that all the sources emit HEPs with 
the same power, it is possible that the power of HECR emission 
depends on sources. In addition, the anisotropic emission of 
HECRs, e.g., HECR generation in relativistic jets, can change 
the contribution of a source to the total flux of HECRs because 
the trajectories of HECRs which can reach the Earth are lim- 
ited. This anisotropy can be also generated by magnetic field in 
the cluster of galaxies even if sources emit HECRs isotropically 
(iDolag et al. 1120091) . 

Heavy nuclei in HECRs are crucial for observing the correla- 
tion between HECRs and their sources due to their larger elec- 
tric charge. The propagation trajectories of heavy nuclei are 
deflected more strongly than those of protons. In this study, we 
focused on the proton component of HECRs. However, the re- 
cent report of the PAO showed that HECRs include a large frac- 
tion of the heavier nucle i based on < X^ax > and the deviations 
of Xmax measurements ("B ellido et al. 1120 0^. while the HiRes 
reported proton-dominated composition IbcIz et al. 2009). On 
the other hand, the PAO simultaneously reported that the cor- 
relation between the arrival directions of HECRs and the po- 
sitions of nearby a strophysical objects dAbraham et al. Il2007l 
l2008HHague et al. II2O O9). In order to consider the consistency 
between the result composition measurement and the correla- 
tion, we briefly investigate the propagation of heavy nuclei in 
OMR 

Heavy nuclei are disintegrated through interactions with 
CMB photons at the highest energies. Since photo- 
disintegration produces lighter nuclei than original nuclei, dif- 
ferent nuclei are contained in detected HECRs even if on ly 
iron nuclei are generated at sources (e.g.. lAllard et al. I (l2005l) ). 
Nevertheless, we consider pure iron composition to see the ef- 
fect of GMF on heavy nuclei clearly. A significant fraction of 
the PAO events has arrival directions directed to Centaurus A 
(Cen A), the nearest radio-loud AGN 3.4 Mpc). Since the 
distance of Cen A is less than the mean free path of photo- 
disintegration of irons at 10^" eV, irons can arrive at the Earth 
from this source without significant photo-disintegration. Also, 
irons can be most easily acceler ated up to 10^" eV following 
the Hillas criterion jHillas Ill984l) . Although pure iron compo- 
sition is an extreme case, a significant fraction of irons can be 
expected if Cen A is really one of the HECR sources. 

Fig. [15] shows examples of the trajectories of irons with the 
energy of 10'^ eV projected onto a plane perpendicular to the 
Galactic disk and including the Galactic center and the Earth. 
Their arrival directions are (£,b) = (180.0°, 0.3°). This figure 
corresponds to the lower right panel of Fig. [T] but the dipole 
component and the turbulent component of GMF are also taken 
into account. The trajectories are complicated and do not follow 
a simple picture explained at Fig. [T] The red and magenta lines 
are trapped by strong field including the turbulent component. 
In any case, the directions of sources are quite different from 
the arrival directions of iro ns. The propaga t ion of iron nuclei in 
GMF is studied in detail bv lGiacinti et al. I ( l2010l) . 




Fig. 15. — Examples of the trajectories of irons with the energy of lO"^ 
eV projected onto a plane perpendicular to the Galactic disk and including the 
Galactic center and the Earth. The position of the Earth is jc = 8.5 kpc. Their 
arrival directions are {E,h) = (180.0°, 0.3°). The BS-S model (red), the BS-A 
model (green), the AS-S model (blue), the AS-A model (magenta), are consid- 
ered as a spiral component of GMF. The dipole component and the turbulent 
component are also taken into account. 

Fig. [T6]demonstrates the arrival directions of the PAO events 
before GMF deflections, assuming that all the PAO events are 
protons (green) or irons (blue). The BS-S model (upper left), 
the BS-A model (upper right), the AS-S model (lower left), the 
AS-A model (lower right) are considered as a spiral component 
of GMF. The dipole component and the turbulent component of 
GMF are also taken into account. These directions correspond 
to the positions of their sources if the effect of IGMF can be 
neglected. The detected arrival directions of the PAO events 
(red) and IRAS galaxies within z = 0.018 (black) are also plot- 
ted. In order to make this figure, the detected PAO events are 
backtracked in the 4 GMF models and their velocity directions 
just outside the Galaxy (at 40 kpc from the Galactic center) are 
plotted. The deflections of irons are quite large, and the corre- 
lation between the arrival directions (before GMF deflections) 
and the galaxy distribution is obviously lost. It is unnatural that 
the deflections by IGMF improve the correlation. Thus, a pure 
iron scenario is problematic in the viewpoint of the spatial cor- 
relation. Since light nuclei including protons may be contained 
in detected HECRs in a real situation, the situations can become 
better, but this figure points out that the contribution of GMF to 
the total deflections of nuclei are not negligible. In this case, 
the fraction of protons (or light nuclei) is an important parame- 
ter. This parameter also gives a hint of cosmic -ray composition 
at sources, i.e., a hint of the generation mechanism of HECRs. 
The arrival directions of irons before GMF modification are de- 
pendent on GMF modelling, as we can see in Fig. [16] Detailed 
modelling of GMF could also improve the correlation between 
the arrival directions of HECRs before GMF modification and 
their sources. 
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Fig. 16. — AiTival directions of the PAO events before GMF deflections. Tlie BS-S model (upper left), BS-A model (upper right), AS-S model (lower left), and 
AS-A model (lower right) are considered as a spiral component of GMF. The dipole component and the turbulent component of GMF are also taken into account. 
The composition of the PAO events is assumed to be protons (green) and irons (bhie). The arrival directions of the PAO events are also plotted (red). Black points 
are IRAS galaxies within 75 Mpc. 
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